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ABSTRACT 

In  studies  of  low-frequency  reverberation  within  the  marine  environment,  a  central  concern  is  the 
relationship  between  reverberation  events  and  morphological  features  of  the  seafloor.  A  time-do¬ 
main  migration  algorithm  for  the  reverberation  intensity  field  is  developed  that  produces  scattering 
coefficient  maps  coregistered  with  a  bathymetry  database.  The  algorithm  is  tailored  to  broadband 
transient  sources  with  good  range  resolution,  and  was  developed  to  analyze  an  extensive  set  of 
reverberation  records  from  a  200-255  Hz  source  collected  on  the  flanks  of  the  Mid- Atlantic  ridge. 
The  precise,  sample-by-sample,  tracking  of  wavefronts  across  elements  of  the  bathymetry  database 
that  forms  the  foundation  of  the  algorithms  implementation  results  in  reverberation  maps  that  show 
a  clear  and  detailed  correlation  between  scattering  and  morphology  with  narrow  scarp  slopes  con¬ 
sistently  highlighted.  Environmentally  induced  asymmetries  in  transmission  loss  and  incidence 
angle  are  exploited  to  break  the  inherent  left-right  ambiguity  of  the  receiver  array.  Iterative  migra¬ 
tion,  assuming  a  dominant  dependence  of  backscatter  on  grazing  angle,  produces  images,  even 
from  individual  records,  that  show  good  ambiguity  resolution.  Results  from  multiple  records  cor¬ 
roborate  the  effectiveness  of  the  ambiguity  resolution  and  demonstrate  the  stability  of  the  scatter¬ 
ing  coefficient  estimates  and  the  acoustic  system. 
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INTRODUCTION 

In  July  of  1993,  a  coordinated  set  of  low  frequency  acoustics  experiments  were  conducted  on 
the  western  flank  of  the  Mid- Atlantic  Ridge  within  a  corridor  that  lay  just  north  of  the  Kane  frac¬ 
ture  zone.  These  experiments  were  designed  to  investigate  and  quantify  the  relationship  between 
acoustic  reverberation  and  seafloor  morphology  and  constituted  the  Main  Acoustics  Experiment 
(MAE)  of  a  Special  Research  Program  (SRP)’  on  bottom  reverberation  sponsored  by  the  Office  of 
Naval  Research.  Included  within  the  experimental  suite  were  monostatic  and  bistatic  reverberation 
measurements  by  the  research  vessels  Cory  Chouest  and  Alliance,  and  near  bottom  measurements 
of  low  grazing  angle  interactions  using  vertical  arrays  of  hydrophones.  One  element  of  the  rever¬ 
beration  experiments  that  makes  them  unique  is  the  availability  of  detailed  bathymetric  and  geo¬ 
physical  data  for  the  area,  collected  on  a  pair  of  SRP  sponsored  large-scale  and  small  scale  geo¬ 
physics  cruises.  The  bathymetric  database  includes  a  map  of  the  entire  SRP  experimental  corridor 
gridded  at  200  m  and  finer  scale  surveys  of  selected  areas  with  resolution  on  the  order  of  10  m  or 
better .  It  is  the  availability  of  these  databases  that  permits  the  investigation,  in  detail,  of  the  corre¬ 
spondence  between  reverberation  returns  and  bathymetry. 

The  process  of  mapping  reverberation  returns  back  onto  the  seafloor  scattering  sites  that  pro¬ 
duced  them  can  be  referred  to  as  charting^  !  Here  we  prefer  to  use  the  term  migration,  because  of 
the  strong  similarity  of  the  methods  employed  here  to  the  Kirchoff  migration  method  used  in  seis¬ 
mic  reflection  processing.  The  principal  difference  between  the  two  is  that  the  reverberation  data  - 
beamformed,  time-series  of  acoustic  intensity  —  are  assumed  to  be  the  result  of  incoherent  scatter¬ 
ing  rather  than  coherent  reflections.  Also,  as  is  the  case  for  seismic  migration,  we  restrict  attention 
to  selected  primary  paths  and  ignore  multiple  reflections  in  the  formulation  of  the  method. 

In  outline,  migration  is  accomplished  in  two  stages;  first  the  seafloor  is  mapped  into  the  time- 
beam  coordinates  of  the  data,  then  the  reverberation  signal  from  a  given  beam  and  small  time 
interval  is  distributed  over  the  corresponding  ensonified  area  of  the  seafloor.  Migration  is  compli¬ 
cated  by  the  fact  that  data  were  recorded  by  a  horizontal  line  array  (HLA)  and  each  beam  is  associ¬ 
ated  with  a  pair  of  directions  that  are  oriented  symmetrically  with  respect  to  the  array  axis,  thus  a 
reverberation  return  could  originate  from  one  of  a  pair  of  distinct  ensonified  areas.  Simple  mapping 
schemes  that  ignore  variations  in  seafloor  depth  produce  reverberation  maps  that  are  perfectly 
symmetric  with  respect  to  the  array  axis  and  provide  no  means  of  resolving,  from  a  single  observa¬ 
tion,  the  inherent  "left-right  ambiguity"  introduced  by  a  HLA.  Accounting  for  the  bathymetry  intro¬ 
duces  "environmental  symmetry  breaking"^’'*  -  bathymetry  induced  variations  in  transmission  loss, 
acoustic  shadowing  -  which  can  be  exploited  by  algorithms,  such  as  the  one  presented  here,  to 
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reduce  substantially  the  left-right  ambiguity,  even  for  a  single  ping.  The  difference  between  the 
environmental  symmetry  breaking  technique  and  the  migration  technique  presented  here  is  that  the 
former  assigns  all  energy  to  either  the  left  or  right  side  if  the  difference  in  transmission  loss  reaches 
a  predefined  threshold,  while  the  current  technique  partitions  energy  according  to  a  set  of  weights 
that  vary  continuously  as  a  function  of  transmission  loss  and  other  factors  such  as  expected  back- 
scatter  strength. 

An  additional  problem  with  simpler  schemes  is  that  the  migrated  images  are  in  the  time-beam 
coordinates  of  individual  transmissions  rather  than  in  global  coordinates  of  the  bathymetry,  com¬ 
plicating  the  comparison  of  multiple  pings  and  studies  of  the  relationship  between  scattering  and 
bathymetry.  The  current  migration  algorithm  avoids  this  problem  by  being  formulated  directly  in 
terms  of  the  backscatter  strength  of  individual  elements  of  the  bathymetry  grid. 

The  algorithm  is  based  on  a  time-domain  formulation  of  the  incoherent  scattering  process  which 
relates  the  scattering  strength  to  the  expected  acoustic  intensity  via  a  large  sparse  matrix.  This 
formulation  means  that  we  could  treat  the  migration  as  a  linear  inverse  problem  and  invert  multiple 
pings  simultaneously  for  scattering  strength.  Combining  pings  with  diverse  look  angles  would 
eliminate  left-right  ambiguity  from  the  solutions.  A  linear  inverse  approach  would  also  allow  regu¬ 
larization  constraints  to  be  employed  on  the  solutions,  an  approach  we  pursue  elsewhere^ 

In  this  paper,  we  solve  the  scattering  equations  approximately  using  a  version  of  iterative 
backprojection  under  the  assumption  that  backscatter  strength  is  spatially  isotropic  and  a  function 
only  of  the  local  ensonification  and  backscatter  angles,  angles  which  are  equal  for  the  monostatic, 
direct  paths  considered  here.  The  advantage  of  using  iterative  backprojection  is  that  it  is 
computationally  fast  and  guaranteed  to  produce  migrated  images  that  satisfy  the  data  exactly  while 
only  weakly  enforcing  assumptions  about  the  behavior  of  scattering  strength.  Iterative  backprojection 
could  potentially  be  used  as  a  basis  for  producing  well  resolved,  migrated  images  in  near  real  time. 
We  use  it  here  to  rapidly  check  the  consistency  of  multiple  images  from  successive  reverberation 
records  of  the  MAE  and  to  produce  preliminary  estimates  of  backscattering  strength  as  a  function 
of  grazing  angle.  The  linear  inverse  approach  outlined  above  will  only  be  effective  if  the  individual 
reverberation  images  are  basically  consistent.  Uncorrected  errors  in,  for  example,  array  heading  or 
source  location  would  substantially  degrade  the  result  of  inverting  multiple  pings. 

We  concentrate  on  the  analysis  of  a  subset  of  the  data,  namely  half  convergence  zone  (1/2  CZ), 
monostatic  reverberations  from  broadband,  200-255  Hz,  linearly  frequency  modulated  (LFM) 
transmissions  recorded  by  the  R/V  Cory  Chouest .  We  chose  this  subset  because  it  has  the  best  range 
resolution  and  simplest  propagation  paths,  and  thus  affords  the  best  opportunity  of  investigating. 
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quantitatively,  the  relationship  between  backscattering  strength  and  bathymetry.  Previous  exami¬ 
nations  of  the  SRP  acoustics  data  have,  for  the  most  part,  concentrated  on  the  longer  range,  continu¬ 
ous  wave  (cw)  data,  and  have  established  that  there  is  a  broad  correspondence  between  high  ampli¬ 
tude  reverberations  and  well  ensonified  features  at  1/2  and  1  1/2  CZ,  with  the  highest  returns  from 
backfacing  ridges.  However,  the  detailed  investigation  and  quantification  of  this  relationship  is  as 
yet  incomplete;  this  paper  represents  one  step  in  that  process. 

I.  OVERVIEW  OF  DATA  AND  EQUIPMENT 

The  MAE  was  conducted  within  an  approximately  4°  by  2°  corridor  on  the  flanks  of  the  Mid- 
Atlantic  ridge  extending  from  45°  W  to  49°  W  and  from  25°  30'  N  to  27°  30’  N.  The  corridor  lies 
within  a  larger  area  designated  as  an  ONR  natural  laboratory®,  which  straddles  the  Mid-Atlantic 
ridge  itself  and  is  bounded  to  the  south  by  the  Kane  fracture  zone.  The  seafloor  fabric  within  the 
corridor  is  dominated  by  lightly  sedimented,  ridge  parallel  abyssal  hills,  but  there  are  also  large 
sedimented  ponds  and  the  area  is  dissected  by  deep  corridors  that  are  the  fossilized,  off-axis  expres¬ 
sion  of  spreading  segment  boundaries  at  the  ridge  axis.  The  feature  whose  backscatter  we  will 
examine  in  detail  here  is  a  ~30  km  long  abyssal  hill,  designated  B’,  that  lies  at  the  western  end  of 
one  of  the  segment  boundary  corridors.  In  particular  we  will  concentrate  on  a  set  of  1/2  CZ,  monostatic 
reverberation  returns  recorded  by  the  R/V  Cory  Chouest  that  span  the  intersection  of  a  pair  of  the 
experimental  runs,  run  5a  and  run  6,  Fig.  1  and  Table  I.  This  intersection  lies  slightly  northwest  of 
a  large  set  of  intersecting  tracks  that  constitutes  the  focus  of  an  extensive  set  of  bistatic  and  monostatic 
experiments  that  had  B’  as  the  primary  target.  We  have  chosen  the  first  intersection  here  because 
the  associated  transmissions  provide  better  ensonification  of  the  front  slope  of  B'  at  1/2  CZ. 

The  basic  features  of  the  R/F  Cory  Chouest  source  and  receiver  arrays  are  shown  in  Fig.  2.  The 
vertical  source  array  consisted  of  10  elements  spaced  at  2.29  m  and  centered  at  101  m,  while  the 
horizontal  receiving  array  consisted  of  128  hydrophone  groups  spaced  at  2.5  m.  Beamforming  of 
the  reverberation  data  was  performed  aboard  ship  using  a  delay  and  sum  beamformer  with  a  Ham¬ 
ming  window,  resulting  in  a  nominal  broadside  beamwidth  of  ~1.3°,  measured  at  the  -3  dB  point. 
The  presence  of  dead  or  improperly  gained  phones  in  the  HLA  resulted  in  nearly  uniform  sidelobe 
level  of  approximately  -30  dB'. 

The  schedule  of  source  transmissions  for  the  experiment  was  conceptually  arranged  in  a  hier¬ 
archical  fashion  with  the  bottom  two  levels  in  the  hierarchy  being  transmission  segments  and  indi¬ 
vidual  source  wavetrains  or  pings.  Segments  were  numbered  consecutively  and  each  one  lasted  12 
minutes.  Ping  numbers  were  also  numbered  consecutively  with  typically  6  source  wavetrains  being 
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transmitted  during  each  segment.  For  the  portions  of  runs  5a  and  6  examined  here,  2  out  of  every  3 
segments  were  used  for  R/V  Coty  Chouest  transmissions  and  the  remaining  segment  for  Alliance 
transmissions.  One  of  the  standard,  usually  the  first,  transmissions  within  a  Cory  Chouest  segment 
was  a  5  s.  duration  linearly  frequency  modulated  (LFM)  signal  chirped  over  the  band  200-255  Hz. 
It  had  the  largest  bandwidth  of  any  of  the  source  wavetrains  and  hence  the  best  range  resolution 
after  matched  filtering,  Ar  =  13.6  m;  it  is  thus  the  best  choice  for  examining,  in  detail,  the  spatial 
structure  of  the  scattering  process. 

The  structure  of  the  LFM  source  wavefield  out  to  beyond  a  1/2  CZ  is  shown  in  Fig.  3.  The 
broadband  transmission  loss  was  calculated  by  integrating  across  the  band,  the  transmission  losses 
predicted  at  individual  frequencies  by  a  wide-angle  parabolic  equation  (PE)  code’.  At  times  during 
the  experiment  the  source  array  was  steered  down  at  various  angles  with  respect  to  the  surface,  but 
for  the  transmissions  of  interest  and  the  calculation,  the  steering  angle  was  set  at  0°.  The  sound 
speed  profile,  displayed  in  Fig.  3a,  was  essentially  the  same  throughout  the  experiment  and  the  sea 
state  was  calm.  The  center  of  the  source  array  at  1 8 1  m  was  in  the  upper  part  of  the  waveguide  with 
the  conjugate  depth  being  at  3800  m  and  the  turning  range,  1/2  CZ  distance,  being  around  33  km. 
The  core  of  the  main  acoustic  beam  contains  a  pair  of  high  amplitude  fringes  which  result  from  the 
interference  between  energy  that  was  up  and  downgoing  at  the  source  level.  Behind  the  core,  the 
two  components  are  sufficiently  separated  in  time,  up  to  ~30  ms,  that  there  is  no  significant  inter¬ 
ference  effect  on  the  amplitude.  Furthermore  in  this  region  the  ray  propagation  angles  diveige  by 
about  2°  (Fig.  4).  Beyond  the  turning  range,  most  of  the  structure  in  the  main  beam  is  due  to 
caustics  in  either  the  up  or  downgoing  source  field.  The  source  field  also  has  a  number  of  sidelobes 
that  intersect  the  seafloor  at  ranges  less  than  15  km. 

li.  MIGRATION  ALGORITHM 

The  goal  of  the  migration  algorithm  is  to  produce  a  map  of  scattering  strength  registered  on  the 
same  grid  as  the  bathymetry  data.  The  basic  observational  data  are  the  intensity  time  series  I((P(,,t) 
associated  with  beam  angle  (p^  that  result  from  squaring  the  matched  filtered  and  beamformed 
pressure  series.  For  direct,  reciprocal  paths  between  source/receiver  and  the  scattering  site,  the 
relationship  between  the  expected  intensity  of  the  backscattered  signal  and  the  scattering  coeffi¬ 
cient  is  taken  as  the  following  integral  over  the  seafloor 

6/t)  =  jj  dAy  g{x,(p  i,(p  (1) 


where 
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Tr,s  are  the  ray  theoretical  travel  times  to  and  from  the  scattering  location.  Y  is  the  product  of  the 
transmission  losses,  'Fg.r,  to  and  from  the  seafloor;  the  gains,  ys,  Yr>  and  Ymf.  associated  with  the 
source  array,  receiver  array  and  matched  filtering;  and  the  source  level,  a  product  of  the  transmis¬ 
sion  level,  Pt,  and  signal  duration,  T.  Since  we  are  dealing  with  broadband  transient  signals,  gains 
for  individual  components  of  the  system  are  defined  in  terms  of  the  ratio  of  integrals  over  the  input 
and  output  pressure  fields*  - 

y-jytpwfyim  (3)  * 

The  scattering  coefficient,  m(x,0,(t))  in  Eq.  (1),  is  a  function  of  position,  x,  and  the  monostatic 
incidence/backscatter  direction  at  the  seafloor  (0,(1)).  The  quantity  g(x,  (pb.  <Pa.  t'Tr  -Ts)  can  be 
termed  the  scattering  function  for  the  system,  it  depends  on  the  details  of  the  scattering  process,  as 
well  as  quadratically  on  the  effective  source  function  at  the  seafloor  and  quadratically  on  the  time 
domain  response  of  the  beamformer  for  a  given  steering  angle  (pb  and  arrival  angle  (pa  at  the  re¬ 
ceiver  array.  The  partition  of  scattering  between  the  scattering  coefficient  and  the  scattering  func¬ 
tion  is  defined  by  normalizing  g(t)  so  that  its  time  integral  is  unity. 

For  direct  paths,  Eq.  (1)  is  a  time  domain  variant  of  the  standard  equation  relating  the  expected 
reverberation  intensity  to  scattering  coefficients  for  cw  signals^*^.  The  additional  element  in  the 
time  domain  formulation  is  the  introduction  of  g,  which  can  be  regarded  as  representing  the  resolu¬ 
tion  function  of  the  acoustic  imaging  system.  Analytic  expressions  or  empirical  forms  could  be 
found  for  g  (and  also  m),  given  explicit  assumptions  about  the  nature  of  the  surface,  for  example 
small  amplitude  roughness^.  In  practice,  we  assume  that  the  details  of  the  time  response  will  not  be 
important  provided  that  the  resulting  scattering  coefficients  estimates  represent  averages  over  a 
sufficiently  large  area.  We  thus  use  a  simple  factored  boxcar  form  for  g 

g((p  i,/9  «/f)  =  bf(p  fl)  =  B{f/f  jB(((p  y-(p  ,)/A(p)  (4) 

where  is  B(t)  is  a  unit  boxcar  centered  at  zero  and  width  one,  and  b((pb,(Pa)  is  the  beamformer 
response,  which  is  in  turn  simplified  to  a  boxcar  of  width  A(p,  the  angular  separation  between 
adjacent  beams.  In  the  limit  tw  =  0,  the  time  dependence,  g*(t),  reduces  to  a  6-function.  A  poten¬ 
tially  more  accurate  representation  of  the  angular  dependence  would  be  to  use  the  response  func¬ 
tion  of  the  beamformer  at  a  representative  frequency  rather  than  a  boxcar . 

We  can  derive  a  sonar  equation  from  Eq.  (1)  by  substituting  for  the  scattering  function  from  Eq. 
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(5) 


(4).  With  g*(t)  set  equal  to  a  S-function,  integrating  from  t,  to  t,+T^  yields 
pX  =  'Pr'Ps  •  Ym/YsYr  pX  m.A 

where 

=  ;{<p^t)  (6) 

and  A  is  the  ensonified  area  corresponding  to  the  integration  limits  in  time  and  beam  angle,  and  it  is 
assumed  that  transmission  losses  etc.  are  constant  over  A.  This  is  equivalent  to  the  following  sonar 
equation  for  the  reverberation  level  R  =  201og  pR 

R  =  -TL, -TL,  +  G  +  M  +  101og{A)  +  S - lOlog  [T,  /  7’)  (7) 

where  TLr,s  are  the  transmission  losses,  G  is  the  combined  system  gains,  M  the  scattering  strength 
and  S  the  source  strength,  all  defined  in  the  obvious  way  from  Eq.  (5).  The  final  term  corrects  for 
differences  between  the  averaging  time  and  the  signal  duration.  Frequently  for  cw  signals  the  aver¬ 
aging  time  is  chosen  equal  to  the  pulse  duration  and  there  is  no  correction.  However,  for  the  LFM 
pulse,  choosing  an  averaging  time  equal  to  the  width  of  the  matched  filtered  pulse,  18  ms,  corre¬ 
sponds  to  a  correction  of  24  dB.  A  problem  with  choosing  a  short  averaging  time  is  that  it  is  then 
hard  to  justify  the  use  of  a  8-function  for  the  scattering  function,  since  effectively  the  convolutional 
nature  of  the  reverberation  process  is  being  ignored. 

The  above  result  is  most  naturally  viewed  as  producing  a  scattering  strength  map  in  the  local 
time-angle  coordinates  of  the  beamformer:  the  scattering  strength  estimate  obtained  from  Eq.  (7)  is 
assigned  to  the  area  A  bounded  by  a  pair  of  isochrons  and  a  pair  of  beam  boundaries.  We  wish, 
instead,  to  produce  a  scattering  strength  map  in  terms  of  a  global  coordinate  system  and  thus  we 
will  use  Eq.  (1)  to  reexpress  the  reverberation  intensity  in  terms  of  contributions  from  individual 
bathymetry  elements.  To  this  end,  we  again  use  the  factored  form  of  the  scattering  function,  Eq.  (4), 
and  assume  that  all  quantities,  except  g*,  associated  with  a  small  patch,  k,  of  the  bathymetric  grid 
are  constant.  The  scattering  response  of  patch  k  can  then  be  approximated  as 

=  Yjt  tnjt  g\t)*  ^  =  Yt  Gkit)  (8) 

where,  for  notational  convenience,  we  have  incorporated  the  beamformer  factor  b((pb,(Pa)  into  the 
general  gain  factor  Yk-  dAk/dt  is  the  rate  at  which  the  patch  area  is  swept  out  by  the  isochrons;  it  is 
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zero  for  t  <  Tmin,k  and  t  >  Tmax,k,  the  minimum  and  maximum  travel  times  associated  with  the 
patch.  Convolution  of  g*  with  the  sweep  rate  function  acts  to  smooth  g*,  suggesting,  once  again, 
that  the  knowledge  of  the  details  of  the  scattering  is  not  critical  provided  that  estimates  are  averages 
over  large  enough  areas,  a  necessary  condition  for  which  is  that  Tmax,k  -  Tniin,k  » tw. 

The  total  scattered  intensity  of  a  single  beam  is  the  sum  of  the  individual  patches  responses,  i.e. 

I{t)  =  |ilk(t)  =  ]SYkmkGk(t)  (9) 

The  summation  implicitly  includes  the  left-right  ambiguity  of  the  receiving  array  through  the  patch 
weighting  Yk,  since  the  beamformer  factor  is  nonzero  only  along  the  conjugate  beam  directions.  As 
before  Eq.  (9)  could  be  explicitly  averaged  by  integration  over  a  time  intervals  T^  which  would 
serve  to  stabilize  the  individual  intensity  values  and  reduce  the  size  of  the  subsequent  matrix  equa¬ 
tions.  If  the  scattering  function  is  taken  as  a  8-function,  averaging  yields 

pl{t)T„  =  Yk  fkit) 

where  is  the  area  of  the  kth  patch  and  f^  is  the  fraction  of  the  area  of  the  kth  patch  ensonified 
during  the  averaging  interval.  If  the  fractional  contributions  are  rounded  to  0  or  1,  then  Eq.  (10)  is 
equivalent  to  Eq.  (5)  of  Ref.  (2).  In  the  latter,  a  synthetic  test  case  was  considered  with  the  source 
wavetrain  taken  as  2s  long  cw  pulse  and  the  averaging  time  taken  equal  to  the  pulse  length,  result¬ 
ing  in  multiple  seafloor  patches  contributing  to  a  single  averaged  intensity  measurement.  In  practi¬ 
cal  applications,  though,  it  would  be  difficult  to  justify  the  use  of  a  8-function  when  the  duration  of 
the  source  pulse  is  long  relative  to  the  size  of  the  individual  patches,  since  implicitly,  as  noted 
above,  it  ignores  the  convolutional  nature  of  the  reverberation  process.  A  scattering  function  whose 
duration  matched  that  of  the  cw  pulse  would  be  more  appropriate,  but  even  so  the  structure  of  the 
function  would  potentially  be  important  for  good  scattering  estimates. 

A  direct  discretization  of  Eq.  (9)  for  time  samples  at  t  =  j  At  is 

Ij  =  ]|lljk  =  iSGjkYk*’^k 

Individual  beam  responses  from  multiple  pings  can  be  combined  to  form  a  single,  large  matrix 
equation  relating  reverberation  to  the  scattering  coefficients  of  the  patches.  For  the  datasets  consid¬ 
ered  below,  the  number  of  scattering  coefficients  would  be  -30,000,  and  the  number  of  intensity 
samples  -20,000,  a  value  which  increases  proportionately  if  multiple  pings  are  combined.  How¬ 
ever,  the  matrix  is  sparse  since  K(j),  the  number  of  patches  illuminated  at  time  sample,  j,  is  on  the 
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of  order  10  for  a  single  ping. 

The  matrix  equation  could  be  solved  as  a  least  squares  problem  with  the  addition  of  regulariza¬ 
tion  constraints  on  the  scattering  coefficients,  using  any  suitable  sparse  matrix  solver  such  as  con¬ 
jugate  gradients*®  or  LSQR**,  an  approach  we  investigate  elsewhere^.  Here  we  choose  to  solve  the 
inverse  problem  approximately  using  a  form  of  iterative  backprojection.  The  advantage  of 
backprojection  is  that  it  is  simple,  fast  and  robust,  while  at  the  same  time  capable  of  providing  a  fair 
degree  of  ambiguity  resolution;  we  use  it  here  to  check  the  ping-to-ping  consistency  of  the  rever¬ 
beration  data.  Although,  backprojection  does  not  have  the  potential  resolution  of  more  complete 
inversion  methods,  especially  for  multiple  pings,  these  methods  will  not  perform  optimally  unless 
the  data  are  essentially  consistent  and  we  can  be  assured  that  there  are  no  unresolved  problems  with 
ship  position  and  array  orientation'*. 

For  a  single  beam,  the  intensity  contribution  assigned  to  a  patch  k  at  time  sample  j  by  the 
backprojection  is  proportional  to  the  weighting  factor  for  that  patch 

Ijk  =  WkIj 

where 


Wk  = 


Yk‘  ^jk 


(13) 


The  backprojection  thus  incorporates  a  degree  of  ambiguity  resolution,  assigning  intensity 
based  on  such  factors  as  transmission  loss  to  the  ensonified  patches  and  also  on  the  basis  of  the  area 
swept  out  -  a  face  pointing  more  directly  into  the  acoustic  beam  will  attract  proportionately  higher 
assigned  intensity. 

The  result  of  the  backprojection  is  an  estimate  of  the  scattering  intensity  function  Ik(t)  associ¬ 
ated  with  each  patch  k.  An  estimate  of  the  scattering  coefficient,  mk,  is  found  by  integrating  the 
scattering  intensity,  Eq.  (9),  and  using  the  fact  that  the  integral  of  g*(t)  is  unity  and  thus  the  integral 
of  G(t)  is  Ak,  the  area  of  patch  k.  The  discrete  time  version  of  the  estimate  is 


mk  = 


^Ijk 

)  ^ 

YkAk 


(14) 


or  substituting  for 
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(15) 


1 

“‘'a  Y.VtT 


This  latter  form  shows  explicitly  that  the  scattering  coefficient  estimate  is,  as  it  should  be,  the  ratio 
of  the  scattered  intensity  to  the  incident  intensity,  scaled  by  the  area. 

The  summation  over  the  individual  intensity  contributions  helps  to  stabilize  the  scattering  coef¬ 
ficient  estimate  for  each  patch.  A  rough  estimate  of  its  expected  variance  can  be  obtained  by  con¬ 
sidering  the  variance  of  the  contributing  intensity  values,  assuming  the  are  approximately  con¬ 
stant.  The  time  bandwidth  product  for  the  LFM  pulse  and  200  m  patches  is  approximately  15,  thus 
the  expected  variance  of  the  sum  is  approximately  1  dB‘^'*. 

If  a  patch  is  illuminated  by  more  than  one  beam,  then  the  scattering  coefficient  estimate  is 
modified  appropriately  to  become  a  weighted  sum  of  the  individual  beam  contributions. 


(16) 


where  is  the  fractional  area  of  patch  k  illuminated  by  beam  n  and  is  the  appropriate  gain 
factor  for  the  beam. 

In  the  examples  examined  later,  we  have  employed  an  iterative  form  of  the  backprojection  in 
which  the  weighting  function  in  Eq.  (12)  is  modified  to 


Wk  = 


Yk- 

^S^)Yk.%.Gk,j 


(17) 


where  ifik  is  the  scattering  coefficient  estimate  from  the  previous  iteration.  The  numerator  is  now 

the  predicted  scattering  intensity  for  each  patch,  and  the  effect  of  the  backprojection,  when  there  is 
a  mismatch  between  the  predicted  and  recorded  intensities,  is  to  modify  each  of  the  assigned  inten¬ 
sities  by  an  equal  decibel  increment. 

The  backprojection  method  has  similarities  to  the  "environmental  symmetry  breaking"  approach 
of  Ref.  3,  particularly  on  the  first  iteration  when  the  factor  controlling  the  assignment  of  energy  is 
usually  transmission  loss.  The  difference  between  the  two  is  that  Ref.  3.requires  the  difference  in 
transmission  loss  to  reach  a  threshold  before  assigning  all  energy  to  one  side,  while  here  assign- 
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ment  of  energy  changes  continuously  as  a  function  of  the  weighting  factors. 

The  above  formulation  is  essentially  a  local  one,  relating  a  deterministic  reverberation  signal, 
consisting  of  sections  of  elevated  intensity,  to  scattering  produced  by  direct  paths  to  and  from  the 
seafloor.  Local  volume  scattering  can  be  considered  to  be  included  as  part  of  the  scattering  signal 
via  the  scattering  function  g.  All  other  contributions  to  the  recorded  reverberations,  including  non¬ 
local,  multiple  scattering  and  long-range,  sub-seafloor  refraction  paths  are  grouped  together  as 
noise.  Forward  modeling  using  a  PE  approach  indicates  that  multipath  contributions  to  the  rever¬ 
beration  field  are  not  significant  for  ranges  below  a  1/2  CZ  except  as  diffuse  arrivals  that  fill  in  the 
noise  floor.  However,  compact  multipath  arrivals  do  become  significant  in  the  time  interval  imme¬ 
diately  beyond  the  direct  path  1/2  CZ  returns,  where  they  fill  what  other  would  be  an  acoustic 
shadow  zone  (13A).  In  principle,  it  would  be  possible  to  incorporate  multiple  scattering  in  the 
above  formulation  by  including  suitable  ray  paths.  However,  such  an  approach  would  quickly  be¬ 
come  unwieldy  if  many  paths  were  included,  and  would  itself  require  an  a-priori  parameterization 
of  the  scattering  process,  which  is  one  of  the  objectives  of  the  study.  Full  wave  solutions,  such  as 
PE,  automatically  handle  the  multipath  problem  but  do  not  circumvent  the  need  for  prior  param¬ 
eterization.  The  validity  of  the  local  scattering  assumption  and  the  possible  influence  of  multipaths 
can  be  examined,  a-posteriori,  by,  determining  whether  there  are  any  regions  of  high  scattering 
coefficient  not  associated  with  identifiable  bathymetric  features  or  that  are  not  stable  between  mi¬ 
gration  images. 

III.  IMPLEMENTATION  DETAILS 

In  implementing  the  migration  algorithm  for  the  SRP  data,  we  used  as  the  grid  for  the  scattering 
coefficient  map  the  200  m  by  200  m  swath  bathymetry  grid  from  the  Hydrosweep  survey.  This  grid 
size  is  smaller  than  the  nominal  cross  track  resolution  of  the  beamformed  data,  which  is  ~  600  m  at 
broadside  at  1/2  CZ,  but  expressed  as  an  ensonification  duration,  270-380  ms,  it  is  considerably 
greater  than  the  18  ms  duration  of  the  compressed  LFM  pulse.  For  such  a  broadband  signal,  the 
compressed  pulse  width  yields  an  optimistic  estimate  of  range  resolution,  since  the  actual  resolu¬ 
tion  would  depend  strongly  on  the  duration  of  the  scattering  response,  including  the  surface  reflec¬ 
tion  delay  at  the  source/receiver.  Whether  the  200  m  grid  is  large  with  respect  to  the  scattering 
function  duration  can  be  checked  a-posteriori  by  assessing  the  stability  of  the  scattering  strength 
maps. 

The  quantities  needed  for  backprojection  of  a  given  ping  are  calculated  in  two  stages,  first 
relevant  propagation  quantities  such  as  travel  time  and  ray  angles  are  found  at  the  vertices  of  the 
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grid,  then  the  sweep  function  G(t),  is  estimated  for  each  200  by  200  m  grid  element,  which  we  refer 
to  here  as  a  patch.  The  spatial  homogeneity  of  the  sound  speed  profile  during  the  MAE,  ignoring 
small  scale  fluctuations,  permits  a  considerable  reduction  in  computational  effort:  propagation  quan¬ 
tities  are  not  found  directly  by  three-dimensional  ray  tracing  but  are  instead  interpolated  from  fields 
calculated  on  regular  two-dimensional,  range-depth,  grids.  This  simplification  to  a  one-dimen¬ 
sional  profile  also  makes  it  computationally  feasible  to  use  PE  calculations,  rather  than  less  accu¬ 
rate  ray  theory,  to  estimate  broadband  transmission  losses.  Two  sets  of  PE  results  were  calculated, 
one  with  a  ten  element  source  array  for  the  forward  direction,  the  other  with  a  single  element  at  the 
HLA  depth  for  the  reverse  direction,  assuming  reciprocity.  The  single  frequency  results  were  then 
integrated  across  the  source  band  to  obtain  the  broadband  TL  estimates  (c.f.  Figure  3).  The  effects 
of  the  free  surface  reflection  are  included  in  the  amplitude  terms  of  the  migration  via  the  TL  calcu¬ 
lation.  However,  only  one  set  of  travel  times  for  the  downgoing  ray  paths  are  calculated,  thus  the 
theoretical  resolution  is  degraded  by  at  least  the  time  difference  in  the  up  and  downgoing  ray  paths, 
~20  ms  (c.f.  Figure  4). 

The  properties  of  the  individual  patches  are  calculated  from  the  vertex  values.  The  patches  are 
assumed  planar,  and  their  orientation  is  found  by  a  least  squares  fit  to  the  vertex  depths.  Similarly 
the  isochron  surfaces  of  two-way  travel  time  are  assumed  to  be  locally  planar  and  are  found  by  a 
least  squares  fit  to  the  travel  times  at  the  vertices.  In  the  following  examples,  we  take  the  scattering 
function,  g(t)  to  be  a  5  function,  and  thus  the  sweep  function,  G(t,)  for  a  patch  is  comprised  of  a 
series  of  straight  line  segments  with  vertices  corresponding  to  times  when  the  wavefront  passes  a 
patch  vertex  (Fig.  5).  From  Fig.  5,  it  is  evident  that  using  a  more  complex  short  duration  function 
rather  than  a  6-function,  for  example  a  1 8  ms  boxcar,  would  have  only  a  minor  effect  on  the  scatter¬ 
ing  coefficient  estimates  after  averaging. 

Geometric  shadowing  calculations  were  performed  prior  to  migration  using  a  mean  bathymet¬ 
ric  profile  along  each  beam  direction.  A  small  transition  region  of  300  m  was  added  to  the  edge  of 
each  ensonified  region  in  order  to  help  prevent  small  numerical  artifacts  or  local  roughness  from 
producing  shadow  zones.  Thus  individual  patches  that  point  away  from  the  local  ray  direction  were 
considered  ensonified,  and  it  was  left  to  the  migration  to  determine  whether  they  produced  small 
backscatter.  During  backprojection,  the  scattering  strengths  of  patches  in  the  shadow  zones  were 
adjusted  at  each  iteration  in  order  to  produce  an  expected  backscattered  intensity  equal  to  the  10th 
percentile  of  the  ensonified  patches.  The  purpose  of  this  was  to  ensure  that  the  preponderance  of 
energy  would  be  assigned  to  an  ensonified  area  if  the  conjugate  side  was  in  shadow,  while  at  the 
same  time  permitting  energy  assignment  if  both  sides  were  in  shadow.  In  this  way  it  was  possible 
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to  keep  track  of  unexpected  scattering  sites  in  the  reverberation  map. 

The  migration  algorithm  was  implemented  primarily  in  MATLAB  but  with  the  sweep  func¬ 
tions  for  individual  patches  calculated  by  an  external  C  function.  The  first  stage  of  the  algorithm, 
the  calculation  of  propagation  quantities  at  the  vertices  of  the  grid  took  6  minutes  of  CPU,  and 
the  second  stage  took  7  minutes  per  iteration  on  a  HP  715/75  workstation  (SPECint95  3.1, 
SPECfp95  3.6).  MATLAB  is  an  interpretted  environment  and  we  conservatively  estimate  that 
these  times  could  be  reduced  by  a  factor  of  four  if  the  algorithm  was  implemented  as  compiled 
code. 

IV.  MIGRATION  EXAMPLE  FROM  B’ 

We  take  as  an  initial  migration  example,  ping  1313  of  segment  415,  for  which  the  RN  Cory 
Chouest  was  located  just  north  and  west  of  the  center  of  the  track  star  near  ridge  B’,  Fig.  1,  and 
had  a  HLA  heading  of  169°,  Table  I.  At  a  1/2  CZ  on  the  starboard  side  is  B’  with  its  backfacing 
scarps  oriented  almost  perpendicular  to  the  beam  directions,  while  at  the  same  range  on  the  port 
side,  a  series  of  N25°E  trending  ridges  point  almost  directly  into  the  beam  directions  (Fig.  1). 
The  reverberation  returns  from  these  ridges  and  the  scarp  of  B’  interfere  to  produce  a  set  of 
crossing,  high  amplitude  events  in  the  data  between  30  and  45  s.  Fig.  6.  When  we  examine  the 
two-way  transmission  loss  to  the  conjugate  areas,  we  find  that  only  the  port  side  ridges  and  the 
upper  part  of  B’  are  shallow  enough  to  project  into  the  center  of  the  acoustics  beam,  although  a 
broad  platform  on  the  port  side  is  also  well  ensonified.  Fig.  7. 

Fig.  8  displays  the  reverberation  map  obtained  after  the  first  iteration  of  the  migration  using 
data  from  20s  out  to  times  corresponding  to  the  1/2  CZ  turning  ranges.  The  reverberation  map  is 
an  intermediate  backprojection  output  and  is  calculated  from  the  sum  of  the  backscattered  en¬ 
ergy  associated  with  each  patch  of  the  grid,  Eq.  (12).  The  reverberation  map  is  a  stable  estimate 
in  the  sense  that  backprojection  is  guaranteed  to  be  intensity  preserving,  and  the  partition  of 
energy  between  simultaneously  ensonified  patches  is  affected  only  by  relative  not  absolute  er¬ 
rors  in,  for  example,  transmission  loss.  From  Fig.  8  we  can  see  that  the  high  amplitude  backscat- 
ter  events  are  predominantly  located  within  areas  that  were  strongly  ensonified.  However,  a 
notable  exception  is  the  strong  scattering  associated  with  a  line  of  small  scarps  on  B'  in  the  region 
(2955-2970, 205-210)  UTM  km,  which  lies  well  outside  the  main  acoustic  beam.  Furthermore, 
it  is  reassuring  that  within  the  strongly  ensonified  regions  the  reverberation  events  tend  to  be 
localized  along  bathymetric  features  with  high  slopes,  such  as  the  scarps  covering  the  upper  parts 
of  B’.  This  impression  is  supported  by  Fig.  9,  which  shows  that  above  a  cutoff  of  about  10°  the 
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mean  scattering  strength  increases  with  local  grazing  angle,  a  dependence  that  arises  even  though 
the  a-priori  assumption  is  that  the  scattering  coefficient  is  constant.  The  cutoff  angle  of  about  10° 
represents  the  point  at  which  the  scattering  strength  curve  reaches  the  noise  floor,  approximately  30 
dB  below  the  maximum  scattering  strength.  The  primary  source  of  noise  is  cross-talk  between  the 
beamformed  channels  and  thus  the  usuable  dynamic  range  of  the  ARSRP  reverberation  data  is  set 
by  the  -30  dB  side-lobe  level  of  the  beamformer.  The  influence  of  the  cross-talk  noise  on  the 
ARSRP  reverberation  data  is  nicely  demonstrated  by  finite-difference  modeling  studies*^. 

Although,  the  initial  migration  produces  a  fair  degree  of  left-right  ambiguity  resolution,  it  is 
clearly  not  perfect  if  local  slope  is  the  sole  or  primary  determinant  of  backscatter  strength.  For 
example,  energy  that  on  the  starboard  side  is  associated  with  the  scarps  of  B',  also  appears  on  the 
port  side  where  it  cuts  across  the  mostly  shallow  slope  bathymetric  fabric.  We  can  test  the  degree  to 
which  slope  accounts  for  strong  backscatter  by  incorporating  the  estimated  angular  dependence  of 
the  scattering  coefficient  into  the  partition  coefficients  for  the  backprojection,  Eq.  (17).  All  events 
associated  with  high  slope  features  will  migrate  solely  to  the  appropriate  left  or  right  image,  only 
anomalous  events  will  remain  ambiguous  and  appear  on  both  images.  For  this  dataset,  the  migra¬ 
tion  effectively  converged  after  three  iterations,  with  only  a  barely  perceptible  change  in  the  mean 
backscatter  coefficient  between  iterations  3  and  4,  Fig.  9.  The  perceived  ambiguity  in  the 
backscattered  energy  images  is  significantly  reduced  and  only  a  few  small  isolated  events  appear 
split  between  the  two  images.  Fig.  10.  The  effectiveness  of  the  ambiguity  resolution  is  more  readily 
assessed  in  the  original  time-beam  coordinates.  Fig.  6.  For  the  most  part  this  view  confirms  the 
impression  that  the  energy  is  cleanly  split  between  the  left  and  right  side.  However,  energy  near  the 
edges  of  some  prominent  reverberation  events  leaks  through  to  the  other  side,  suggesting  that  there 
may  be  a  slight  error  in  the  array  orientation  or  that  the  migrated  image  could  be  improved  by 
incorporating  some  cross  talk  between  beams  in  the  beamformer  factor  b,  Eq.  (4). 

The  scattering  strength  map  for  the  4th  iteration  is  displayed  in  Fig.  11.  As  expected  the  elon¬ 
gated  scarp  slopes  on  B’  are  highlighted  in  the  image,  but  compared  to  the  energy  image,  the 
prominence  of  the  port  ridges  is  reduced,  their  previous  prominence  being  due  solely  to  strong 
ensonification.  As  was  to  be  expected  from  the  mean  behavior,  there  is  an  overall  correspondence 
between  the  grazing  angle  and  scattering  strength  with  the  normalization  emphasizing  the  scatter¬ 
ing  strength  of  the  lower  slopes.  Also,  predictably,  the  scattering  strength  map  is  noisier  than  the 
reverberation  one,  since  now  global  as  well  as  local  errors  in  estimated  quantities  are  important. 
For  example,  any  mismatch  between  the  predicted  and  actual  pattern  of  transmission  loss  will  be 
reflected  in  the  scattering  strength  map. 

Final  Report  ONR  Grant  N00014-96-1-1004  "Numerical  Studies  of  Seafloor  and  Subseafloor  reverberation”  14 


V.  COMPOSITE  IMAGES  OF  B’ 

We  can  assess  the  validity  of  the  single  ping  ambiguity  resolution  and  simultaneously  gain 
insight  into  the  general  quality  of  the  data  by  comparing  the  migration  results  from  multiple 
ensonifications  of  the  same  target  with  nearly  the  same  ship  location  but  different  headings.  The 
standard  means  of  obtaining  multiple  looks  at  a  target  is  to  use  reverberation  records  from  crossing 
paths  and  this  we  do  by  comparing  images  from  runs  5a  and  6  of  the  MAE  which  cross  at  an  angle 
of  approximately  25°.  But  we  can  also  gain  useful  information  from  multiple  pings  from  a  single 
run,  taking  advantage  of  small  changes  in  the  heading  of  the  receiver  array  between  pings.  If  a 
reverberation  event  is  migrated  to  its  true  origin  on  the  seafloor,  then  its  location  will  not  change 
when  the  array  heading  changes  by  0.  Conversely  a  false  scatterer  that  resulted  from  migrating 
energy  to  the  wrong  side  of  the  array  will  rotate  through  20  since  it  is  located  at  the  mirror  image  of 
true  scatterer  in  the  array.  Fig.  12.  Thus  the  incorrectly  migrated  energy  will  move  in  a  predictable 
fashion  in  step  with  the  changes  in  array  orientation.  This  movement  can  most  readily  be  appreci¬ 
ated  by  combining  the  migration  results  from  successive  pings  into  an  animation  sequence.  For  the 
5  consecutive  LFM  pings  from  run  5a  considered  here,  the  maximum  difference  in  array  heading  is 
3°  which  is  sufficient  to  move  false  scatterers  across  3-4  beams  in  the  migrated  image. 

We  have  investigated  the  consistency  of  the  migrated  images  using  5  consecutive  broadband 
LFM  pings  from  both  run  5a  and  run  6  spanning  their  intersection  at  26°  36'  N,  -47°  47'  W  ( (2946, 
222),  Fig.  1 .).  Fig.  1 3  displays  the  scattering  strength  maps  for  the  northern  foreslopes  of  B  ’  for  9  of 
the  pings.  Perhaps  the  most  striking  aspect  of  the  images  is  that  relatively  narrow  scarp  slopes 
appear  consistently,  with  only  relatively  minor  ping-to-ping  variation  in  strength  and  location.  The 
images  from  run  6  (pings  1683-1696)  are  cleaner  than  the  ones  from  run  5a  (pings  1313-1336)  due 
to  the  fact  that  the  conjugate  area  on  run  6  is  the  source  of  much  less  backscatter.  For  run  5a,  the 
lower  portion  of  the  images  below  2950  includes  residual  mislocated  energy  from  the  port  side 
ridges.  A  port  side  origin  is  supported  by  the  fact  that  it  does  not  appear  in  the  images  from  run  6 
and  that  its  location  moves  between  the  successive  run  5a  images.  Similarly  energy  that  appears  in 
the  run  5a  images  near  the  center  of  a  small  basin  at  (2960,202)  can  also  be  attributed  to  a  partial 
failing  of  the  ambiguity  resolution. 

Averaging  the  individual  images  effectively  suppresses  false  scatterers  attributable  to  energy 
leakage  and  multipaths  and  emphasizes  the  consistency  of  direct  path,  high  backscatter  from  the 
narrow  ridges.  Fig.  14.  All  the  high  backscatter  regions  correspond  to  scarp  features  that  are  ensonified 
at  high  local  grazing  angles.  However,  the  converse  is  not  true,  there  are  features,  most  notably  the 
front  edge  of  a  small  platform  centered  (2955, 198)  that  although  illuminated  at  high  angle  does  not 
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produce  high  backscatter.  Spots  of  high  backscatter  do  appear  along  the  front  edge  of  the  platform 
in  the  individual  migration  images,  perhaps  indicative  of  backscatter  speckle  that  fluctuates  in 
response  to  changes  in  the  ensonification.  However,  since  the  strongest  backscatter  events  appear 
in  the  noisier  run  5a  images  it  is  more  likely  that  the  high  ensonification  angles  succeed  in  captur¬ 
ing  port  side  energy  during  migration. 

For  each  of  the  10  pings,  we  have  used  the  migration  results  to  estimate  the  variation  of  mean 
backscatter  strength  with  grazing  angle.  Fig.  15.  The  backscatter  curve  was  found  using  the  Loess 
algorithm*^  assuming  that  the  local  grazing  angles  were  correct.  The  Loess  algorithm  finds  a  point 
on  the  backscatter  curve  for  a  given  grazing  angle  by  fitting  a  line  to  a  subset  of  the  data  centered 
on  the  grazing  angle  using  weighted  least  squares.  Formal  confidence  intervals  for  the  mean  were 
estimated  from  200  bootstrap  samples'^  for  each  ping,  and  were  typically  1-2  dB  at  the  95  % 
confidence  level,  increasing  at  high  grazing  angles.  Over  the  interval  from  10°  to  45°,  the  mean 
scattering  strength  curves  for  run  5a  are  roughly  comparable  to  a  Lambert’s  law  type  behavior  with 
a  Mackenzie  parameter  of  -27  dB,  although  the  curves  are  distinct  based  on  the  formal  confidence 
intervals.  The  backscatter  estimates  for  a  given  grazing  angle  are  distributed  approximately  nor¬ 
mally  with  a  standard  deviation  of  4-5  dB,  Fig.  16. 

The  split  is  probably  a  consequence  of  differences  in  the  distribution  of  the  relatively  small 
number  of  patches  that  are  ensonified  at  high  grazing  angles.  For  run  5a,  73%  of  the  700-900 
patches  ensonified  at  angles  greater  than  25°  (approximately  3%  of  the  total)  come  from  the  B’ 
scarps,  while  for  run  6,  this  percentage  drops  to  63%  with  a  greater  number  of  high  angle  targets  on 
the  starboard  side  that  are  small  and  relatively  isolated.  When  the  backscatter  curves  are  split  into 
port  and  starboard  contributions,  the  gap  between  the  starboard  side  curves  that  includes  B’  is 
reduced,  although  not  eliminated.  Moreover,  the  port  side  targets  display  a  weaker  dependence  on 
angle.  Fig.  17,  and  are  consistent  over  the  two  runs.  The  difference  in  scattering  strength  may 
indeed  reflect  an  intrinsic  difference  in  the  backscatter  level  between  the  areas.  However,  it  may 
also  be  a  consequence  of  processing  artifacts  and  of  the  limitations  of  the  bathymetry  database.  The 
size  of  the  port  side  targets  tend  to  be  near  or  below  the  angular  resolution  of  the  system,  and  thus 
the  intrinsic  averaging  of  the  beamforming  coupled  with  any  small  errors  in  the  heading  calcula¬ 
tions  will  tend  to  mute  the  influence  of  high  angle  patches.  Also  the  number  of  high  angle  targets  is 
small  enough  that  imperfections  in  the  bathymetry  database  such  as  small  noise  spikes  that  produce 
erroneous  high  angle  targets  would  bias  the  scattering  strength  estimates  downward. 
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VI.  CONCLUSIONS 

We  have  developed  an  algorithm  for  migrating  reverberation  data  that  is  based  upon  time- 
domain  expressions  for  the  expected  intensity  of  the  reverberation  field.  The  time  dependent  as¬ 
pects  of  the  field  are  parameterized  by  a  function  referred  to  as  the  scattering  function,  which 
includes  both  stochastic  components  dependent  on  the  details  of  the  seafloor  scattering  and  deter¬ 
ministic  ones  due  to  the  characteristics  of  the  acoustic  imaging  system.  We  argue  that  the  details  of 
the  scattering  function  are  not  critical  provided  that  scattering  coefficients  are  estimated  over  patches 
of  the  seafloor  with  response  times  that  are  long  compared  with  the  characteristic  time  of  the  scat¬ 
tering  function,  or  equivalently  whose  area  is  large  compared  to  the  correlation  scales  of  the  scat¬ 
tering  process®.  This  approach  is  comparable  to  reducing  speckle  in  optics  by  averaging  over  a 
scanning  aperture’®.  Time  averaging  has  the  additional  advantage  of  producing  stable  estimates 
from  individual  reverberation  records,  the  trade-off  is  that  spatial  resolution  is  reduced.  However, 
even  with  averaging,  the  broadband  LFM  pulse  of  the  SRP  experiments  is  capable  of  yielding 
consistent  scattering  strength  estimates  at  the  resolution  of  the  main  bathymetry  database, 

A  key  feature  of  the  algorithm  is  that  the  results  are  computed  in  the  global  coordinate  system 
of  the  bathymetry  database  rather  than  the  beam-travel  time  coordinates  of  individual  pings,  facili¬ 
tating  the  comparison  of  results  from  multiple  pings.  The  basis  for  the  coordinate  conversion  is  the 
precise,  sample-by-sample,  tracking  of  the  acoustic  wavefronts  across  the  individual  patches  and 
the  assumption  that  the  scattering  coefficient  is  constant  within  individual  elements,  patches,  of  the 
bathymetry  database.  These  two  assumptions  allow  us  to  reexpress  the  reverberation  response  as  a 
large  sparse  matrix  equation  with  the  scattering  coefficients  as  the  unknowns.  The  matrix  formula¬ 
tion  enables  multiple  records  to  be  combined  simply  into  a  single  linear  inversion  problem. 

In  this  paper  we  have  taken  a  processing  oriented  approach  taking  advantage  of  the  sparsity  to 
solve  the  matrix  equation  approximately  using  iterative  backprojection  rather  than  attempting  to 
solve  the  equation  directly  as  part  of  a  linear  inversion.  We  describe  this  process  as  migration  since 
the  method  is  analogous  to  the  Kirchoff  migration  procedure  of  seismic  processing,  the  principal 
difference  being  that  here  it  is  applied  to  incoherent  scattering  rather  than  coherent  reflections.  It 
has  the  processing  virtue  of  producing  reverberation  and  scattering  coefficient  maps  rapidly  thus 
facilitating  efficient  analysis  of  multiple  records.  The  reverberation  maps,  which  are  produced  by 
backprojection  as  an  intermediate  output,  represent  an  exact  partition  of  the  record  intensity  and  are 
thus  guaranteed  to  be  stable  and  allow  us  to  assess  directly  the  effectiveness  of  the  left-right  ambi¬ 
guity  resolution. 

The  performance  of  the  migration  algorithm  is  illustrated  using  two  sets  of  broadband  monostatic 
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reverberation  records  from  a  pair  of  crossing  runs  of  the  SRP  experiment  on  the  flanks  of  the  Mid- 
Atlantic  ridge.  At  a  1/2  CZ,  iterative  migration  produces  well  resolved  maps  of  reverberation  and 
scattering  strength.  In  the  reverberation  maps,  there  is  a  clear  correspondence  between  the  location 
of  high  amplitude  reverberation  events  and  bathymetric  features  with  either  low  transmission  loss, 
high  local  grazing  angle  or  both.  In  particular,  narrow  scarps  less  than  1  km  wide  are  highlighted 
consistently.  The  migration  also  produces  a  plausible  resolution  of  the  left-right  ambiguity,  divid¬ 
ing  the  majority  of  events  in  a  reverberation  record  cleanly  between  the  two  ensonified  areas.  The 
migration  algorithm  can  be  regarded  as  an  environmental  symmetry  breaking  method  that  uses 
seafloor  induced  asymmetries  to  resolve  the  inherent  left-right  ambiguity  of  the  receiving  array. 
Although  the  migration  images  are  by  no  means  unique,  the  improvement  of  the  event  separation 
with  iteration  and  the  clear  correspondence  between  events  and  bathymetric  features  provides  strong 
a-posteriori  justification  for  the  inclusion  of  mean  backscatter  strength  as  a  factor  in  the  backprojection 
weights. 

The  results  of  multiple  migrations  from  a  single  run  and  from  crossing  runs  corroborates  the 
overall  effectiveness  of  the  left-right  ambiguity  resolution  for  individual  records,  and  supports  the 
idea  that  scattering  strength  is  primarily  a  function  of  grazing  angle.  The  majority  of  high  scattering 
strengths  areas,  primarily  associated  with  scarp  slopes,  appear  consistently  in  all  migration  images. 
A  few  potentially  anomalous  areas  move  in  step  with  changes  in  array  heading  within  a  single  run 
and  tend  to  be  absent  in  images  from  the  crossing  run.  These  areas  are  identified  as  energy  misassigned 
to  the  wrong  side  of  the  array,  and  are  attenuated  when  multiple  maps  are  stacked  together.  Em¬ 
ploying  multiple  heading  directions  from  crossing  runs  is  the  standard  means  of  reducing  left-right 
ambiguity.  The  present  results  demonstrate  that  even  small  heading  changes  (e.g  2-3°)  from  a  single 
run  can  be  profitably  employed  to  reduce  ambiguity.  The  consistency  of  the  migration  images 
demonstrates  the  stability  of  the  acoustic  imaging  system  used  for  the  SRP  experiments  and  is  a 
good  indication  that  linear  inversion  could  be  applied  successfully  to  the  data. 

The  migration  results  indicate  that  transmission  loss  and  grazing  angle  arethe  dominant  factors 
controlling  reverberation  response  in  the  study  area.A  dependence,  in  the  mean,  of  backscatter 
strength  on  grazing  angle  appearsafter  the  first  iteration  ofthe  migration  when  the  backprojection 
weights  contain  no  bias  towardspatches  with  high  grazing  angles.  We  thus  conclude,  as  have 
previousinvestigators^’^’^”^,  that  the  correlation,  in  the  mean,  between  backscatter  strength 
andgrazing  angle  is  a  robust  feature  of  the  data. 

Weaker  dependencies  of  scattering  strength  on  spatial  and  azimuthalvariations  in  scattering 
strength  were  not  explicitly  sought  for  in  thepresent  analysis.  However,  the  difference  in  the  mean 
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scattering  curvesbetween  B’  and  the  conjugate  sides,  Fig.  17, indicate  that  there  may  be  an  azi¬ 
muthal  dependence  resolvable  in  the  data.The  prominent  ridges  of  B’  are  ensonified  almost  perpen¬ 
dicular  to  thestructural  grain,  whereas  the  starboard  side  ridges  with  the  weakerbackscatter  are 
ensonified  parallel  to  thegrain.  Resolution  of  questions  such  as  azimuthal  dependence 
ofbackscattering  or  the  influence  of  sediments  on  backscattering  will  requirefurther,  more  detailed 
analysis  of  the  data . 

Above  a  noise  floor  of  about  5°,  the  estimated  dependence  of  the  mean  scattering  strength  on 
grazing  angles  is  approximately  equal  to,  but  formally  statistically  distinct  from,  that  of  Lambert’s 
law  with  a  Mackenzie  parameter  of  -27  dB.  Although  the  general  magnitude  of  the  increase  in 
backscatter  strength  appears  robust,  the  detailed  dependence  on  grazing  angle  must  be  treated  with 
some  caution  as  the  accuracy  of  the  slopes  derived  from  the  swath  bathymetry  is  itself  questionable, 
especially  when  the  high  slopes  are  associated  with  smaller  features  near  the  resolution  limit  of  the 
system'®. 
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segment 


Location  (Deg)  Location 

Lat  Lon  Lat 


HLA  Array 
Heading 


a 

415 

1313 

26°  38.6’  N 

-47°  48.6’ W 

a 

417 

1318 

26°  37.4'  N 

-47°  48.1’ W 

2947.8 

221.1 

a 

420 

1327 

26°  35.7'  N 

-47°  47.5’ W 

2944.7 

221.9 

a 

421 

1331 

26°  35.2'  N 

-47°  47.3’ W 

2943.7 

222.3 

a 

423 

1336 

26°  34.0'  N 

-47°  46.9’ W 

2941.4 

222.9 

547 

1683 

26°  38.7’  N 

-47°  47.3’ W 

2950.2 

222.5 

548 

1687 

26°  38.1’ N 

.470  47.4’ W 

2949.1 

222.2 

) 

550 

1692 

26°  36.9'  N 

-47°  47.6’ W 

2946.9 

221.8 

) 

551 

1696 

26°  36.3’  N 

.47°  47.7’ W 

2945.7 

221.6 

) 

554 

1705 

26°  34.5'  N 

-47°  48.1’ W 

2942.4 

220.9 

X  v/u.vr 
166.1 

169.5 
169.1 

167.6 

197.3 
195.9 
194.8 

194.6 

195.3 


Table  I:  List  of  reverberation  records,  pings,  analyzed  in  this  paper 
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Figure  Captions 


FIG.  1.  Bathymetry  and  a  portion  of  R/V  Cory  Chouest  track  lines  (heavy  black  lines)  in  the  vicin¬ 
ity  of  target  B',  an  elongated  ridge  at  the  western  edge  of  the  Fig.  In  this  paper  we  analyze  rever¬ 
beration  data  from  Pings  at  the  start  of  run  5a  and  6  of  the  MAE  (heavy  white  circles). 

FIG.  2.  Sketch  of  the  source  and  receiver  array  of  the  R/V  Cory  Chouest.  The  10  element  source 
array  had  a  spacing  of  2.29  m  while  the  126  element  receiver  array  had  a  spacing  of  2.5  m. 

FIG.  3.  (a)  Representative  up  (blue)  and  downgoing  (red)  rays  paths  within  the  main  acoustics 
beam  from  a  point  source  at  181  m,  the  center  depth  of  the  vertical  source  array.  The  reference 
sound  speed  (black)  has  a  surface  velocity  of  1540  m/s  and  a  minimum  velocity  of  1494  m/s  at  1.1 
km.  The  downgoing  ray  paths  are  simple  but  caustics  form  beyond  the  1/2  CZ  at  30-35  km.  (b)  PE 
calculation  of  the  tranmission  loss  for  a  broadband,  200-255  Hz,  source  from  the  10  element  source 
array.  Transmission  loss  variations  near  the  axis  of  the  downgoing  main  beam  are  the  result  of 
interference  of  the  downgoing  and  surface  refracted  energy.  Transmission  loss  is  more  variable 
beyond  a  1/2  CZ  due  to  the  presence  of  caustics. 

FIG.  4.  (a)  Grazing  angle  of  rays  at  a  depth  of  3.5  km.  Arrivals  are  for  rays  from  the  top,  middle  and 
bottom  elements  of  the  source  array.  There  is  approximately  a  2°  difference  in  ray  angle  between 
surface  reflected/refracted  rays,  upper  cluster,  and  rays  leaving  the  source  travelling  downwards. 
At  this  depth,  the  penumbra  of  the  main  beam  starts  at  25  km  and  extends  to  about  28  km.  (b) 
Difference  in  travel  times  at  3.5  km  for  up  and  downgoing  rays  at  the  source. 


FIG.  5.  Example  sweep  area  function  for  a  seafloor  patch.  The  sweep  function  consists  of  linear 
segments  with  joints  corresponding  to  the  passage  of  the  beam  wavefront  (dash  lines)  past  the  patch 
vertices.  For  the  calculation  both  the  patch  and  wavefront  are  assumed  locally  planar.  (Patch  area  is 
normalized  to  unit  area). 
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FIG.  6.  (a)  Reverberation  data  for  ping  1313,  broadband  LFM  transmission  of  s415.  The  displayed 
data  spans  the  time  interval  of  reverberation  returns  from  1/2  CZ  and  beam  numbers  15-120,  which 
avoids  endfire  (Beam  nos.  64/65  are  the  broadside  beams).  To  facilitate  comparision  with  the  rever¬ 
beration  maps,  the  data  has  been  averaged  over  an  interval  of  270  ms  and  corrected  for  the  average 
number  of  ensonified  patches,  a  correction  of  -11.1  dB.  The  prominent  returns  between  30  &  45s 
are  a  combination  of  returns  from  scarps  on  the  front  face  of  B'  on  the  starboard  side  of  the  array 
and  ridges  on  the  port  side  that  point  almost  directly  into  the  beams,  (b)  &  (c)  estimated  port  and 
starboard  time  series  after  iterative  migration.  The  split  is  plausible  and  for  the  most  part  clean, 
although  some  residual  energy  appears  to  be  misplaced. 

FIG.  7.  Two  way  transmission  loss  for  ping  1313  projected  onto  the  port  and  starboard  side  bathym¬ 
etry.  Strongly  ensonified  features  include  the  top  of  B'  on  the  starboard  side,  and  the  high  ridges 
pointing  approximately  into  the  beams  on  the  port  side.  The  uppermost  parts  of  the  two  prominent 
ridges  on  the  port  and  starboard  side  colored  white  project  up  above  the  acoustic  beam  and  are  thus 
not  included  in  the  calculation.  The  current  figure  does  not  take  into  account  geometric  shadowing. 

FIG.  8.  Reverberation  map  from  the  initial  migration  of  ping  1313.  The  left-right  ambiguity  reso¬ 
lution  is  based  primarily  on  ensonification  level,  although  some  dependence  on  local  incidence 
angle  is  included  via  the  sweep  rate  function.  High  amplitude  returns  due  to  local  slope  such  as 
those  from  the  scarp  slopes  of  B'  are  not  fully  resolved  and  also  appear  in  the  conjugate  image. 

FIG.  9.  Estimate  of  backscatter  strength  as  a  function  of  grazing  angle  for  the  4  migration  iterations 
of  Ping  1313.  Iteration  1-  dotted,  2-  dash-dot,  3-dashed,  4-solid.  For  iteration  1  there  is  no-apriori 
assumption  of  increased  backscatter  with  grazing  angle,  yet  the  mean  backscatter  strength  shows  a 
correlation  with  grazing  angle  down  to  the  noise  floor  at  about  10°.  For  subsequent  iterations  the 
mean  backscatter  of  the  previous  iteration  is  used  and  convergence  occurs  in  3-4  iterations. 

FIG.  10.  Reverberation  map  of  ping  1313  after  four  iterations  .  Compared  to  iteration  1,  the  left- 
right  ambiguity  resolution  is  significantly  improved,  with  a  much  weaker  ghosting  of  the  B'  scarps 
in  the  port  side  image. 

FIG.  ll.Scattering  strength  maps  derived  from  the  reverbearation  map.  Fig  10.,  by  correcting  for 

transmission  loss,  area,  etc.  The  previous  prominence  of  the  top  of  B'  and  the  port  side  ridges  is 

revealed  as  primarily  a  consequence  of  high  ensonification  levels,  while  the  B'  scarp  slopes  are 

identified  clearly  as  regions  of  high  backscatter. 
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FIG.  12.  With  incomplete  ambiguity  resolution,  energy  from  a  scatterer  will  also  be  migrated  to  a 
false  scatterer  location;  that  is,  the  mirror  image  in  the  line  array  of  the  true  location.  These  loca¬ 
tions  will  remain  paired  for  different  pings  (black  circles)  if  the  ship  heading  remains  constant,  but 
the  location  of  the  false  scatterer  will  rotate  through  20  when  the  array  heading  changes  by  0. 


Fig.  13.  Scattering  coefficient  maps  of  a  subarea  of  the  B’  abyssal  hill  for  5  consecutive  broadband 
LFM  pings  of  run  5a,  1313  through  1336,  and  4  consecutive  pings  for  run  6, 1683-1696.  Most  of 
the  scarp  slopes  in  the  area  show  up  consistently  as  regions  of  high  backscatter.  A  notable  exception 
is  the  edge  of  a  small  platform  near  (2955,200)  which  is  only  intermittently  highlighted  possibly 
only  as  the  result  of  being  mistakenly  assigned  energy  from  the  conjugate  area. 


Fig.  14.  Comparison  of  the  mean  backscatter  strength,  left  panel,  found  by  stacking  the  10  migrated 
pings,  with  the  mean  illumination  angle,  right  panel,  expressed  in  terms  of  the  cosine  of  the  local 
grazing  angle.  For  the  most  part,  there  is  a  close  correspodance  between  high  backscatter  targets 
and  local  grazing  angle:  all  prominent  backscatter  targets  correspond  to  scarp  features  ensonified  at 
high  grazing  angles.  However,  there  are  some  features,  most  notably  the  front  slope  of  a  small 
platform  centered  (2955,  198)  that  although  illuminated  at  high  angle  does  not  produce  a  high 
backscatter. 

FIG.  15.  Estimate  of  scattering  strength  vs.  grazing  angle  for  pings  from  run  5a  (solid)  and  run  6 
(dashed)  with  bootstrap  estimates  of  +2s  errors  in  the  mean.  For  reference,  Lamberts  law  curve 
using  a  Mackenzie  parameter  of  -27  dB  is  included,  dotted  curve. 

FiG.  16.  Representative  histogram  of  scattering  strength  relavite  to  mean  trend  from  ping  1313  for 
a  bin  extending  from  21°-29°.  There  are  1385  patches  in  the  bin.  The  distribution  is  approximately 
gaussian  with  a  =  4.4  dB 

FIG.  17.  Port  and  starboard  estimates  of  mean  scattering  strength  for  thelO  pings.  Estimates  from 
run  5a  are  solid  lines  and  those  from  run  6  are  dashed.  The  estimates  from  the  starboard  side  that 
ensonify  B’  show  a  consistently  stronger  dependence  on  grazing  angle  and  there  is  less  of  a  gap 
between  run  5a  &  6  estimates. 
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Figure  1 
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Figure  2 
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Figure  3 
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Figure  5 
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Figure  6 
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Figure  7 
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Figure  10 
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Figure  12 
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%  Matlab  tn-file  to  take  the  derivative  of  a  geoid  height 
%  profile  and  compare  a  first  difference  derivative  to  the 
%  FFT  derivative  approach. 

% 

%  The  file  geosatd.dat  is  loaded. 

%  The  first  column  of  the  file  is  time  (not  used  here) . 

%  The  second  column  of  the  file  is  latitude  in  degrees. 

%  The  third  column  of  the  file  is  longitude  in  degrees. 

%  The  fourth  column  is  geoid  height  in  meters. 

%  The  fifth  column  is  gravity  anomaly  in  milligals. 

%  The  sixth  column  is  uncertainty  in  gravity  anomaly. 

%  (note  this  is  a  descending  track  so  latitude  decreases 
%  with  increasing  time.) 

% 

%  There  are  512  evenly  spaced  points.  The  total  length  of  the 
%  profile  is  1705000  m. 

% 

%  Get  the  hwlpart.m  file  and  the  data  file  using  fetch  or  ftp 
%  baltica.ucsd.edu 
%  cd  pub/class/hwl 
%  mget  * 

% 

load  geosatd.dat 
lat=geosatd ( : , 2 ) ; 
lon=geosatd ( : , 3 )  ; 
geoid=geosatd ( : , 4 ) ; 

% 

%  remove  the  mean  and  window  the  profile  to  minimize  edge  effects 
% 

nx=length(lat) ; 

window=hanning(nx) ; 

geoid= (geoid-mean (geoid) ) . *window; 

% 

%  plot  geoid  height 
% 

subplot (2,1,1),  plot (lat, geoid) ; 
ylabel (' geoid  height  (m)  ') 

% 

%  get  the  profile  length  and  data  spacing 
% 

L=1705000. ; 
dx=L/nx; 

% 

%  compute  the  derivative  using  the  first  difference  formula. 

%  multiply  the  dimensionless  slope  by  le06  to  get  microradians. 

%  duplicate  the  last  point  to  extend  the  length  from  511  to  512. 

%  note  there  is  a  shift 
% 

vd  =  1 . e6*diff (geoid) /dx; 
nd  =  length (vd) ; 
vd  =  [vd ;  vd  (nd :  nd)  ]  ; 

% 

%  now  compute  the  derivative  using  the  derivative  property. 

%  first  compute  FFT  and  do  an  fftshift  to  make  it  easy  to  gerenate 
%  wavenumbers . 

% 

eg  =  f ftshift (f ft (geoid) ) ; 

% 

%  generate  the  wavenumber.  if  these  wavenumbers  are  wrong  the 


